rm(list=ls(all=TRUE))  # new
library(pscl)

subset = 0 # 0: all data, 1: small subset,  101: 101st congress


if (0) {  # 1: load and create dta,  0: load saved one  
	##loading data
	if (0) {  # 1: load csv, 0: load Rdata 
		dta1 = read.csv("Rollcallmatrix101to109corrected.csv", na.strings = c("NA","#N/A"))
		dta2 = read.csv("Rollcallmatrix110to112corrected.csv", na.strings = c("NA","#N/A"))
		dta1 = dta1[,-ncol(dta1)]  # the last two columns are empty
		dta1 = dta1[,-ncol(dta1)]  # the last two columns are empty
		dta2 = dta2[,-ncol(dta2)]  # the last column is empty
		z = which(dta1[,1]!=dta2[,1])
		if (length(z)>0) {
			print(z)
			print("row mismatch")
		}
		dta = cbind(dta1, dta2[,-1])  # the first column of dta2 is the same as that of dta1
		save(dta, file="raw_dta.Rdata")
	} else {
		load("raw_dta.Rdata")
	}

	## remove unanimous votes
	drop.unani <- function(){
		z = matrix(FALSE, ncol = 2, nrow = ncol(dta))
		for (i in 1:ncol(dta)){
			z[i,] = c(all(dta[,i]==0,na.rm = TRUE), all(dta[,i]==1, na.rm=TRUE))
			if (all(is.na(dta[,i]))) { cat("Column ",i,"(",colnames(dta)[i],") is all NA.\n") }
		}
		unani = (z[,1] | z[,2])
		cat("total",dim(dta)[2],"votes\n")
		cat(sum(z[,2]),"(yay) &",sum(z[,1]),"(Nay) unanimous votes\n")
		cat(sum(unani),"votes dropped\n")
		dta[,!unani]
	}

	check.coding <- function(dta){
		for (i in 1:ncol(dta)){
			z=which(!(dta[,i]==1 | dta[,i]==0 | is.na(dta[,i])))
			if (length(z)>0) {
				cat("column ",i,"\n")
				print(z)
			}
		}
	}

	members =as.vector(dta[,1])
	dta=dta[,c(-1)]
	check.coding(dta)
	dta=drop.unani()

	dta=as.matrix(dta)

	data.list = list(dta=dta, members=members)

	save(data.list, file = "dta.Rdata")

} else {
	load("dta.Rdata")
	dta = data.list$dta
	members = data.list$members
}

#creating subset of dta if necessary
if (subset==1) {
	#creating subset of dta
	dta = dta[1:20,1:100]

	members = members[1:20]
} else if (subset==101) { # (101th Congress)
	dta = dta[1:548,1:1353]   # colnames(dta)[1353:1360] #votes in the 101th congress
	members = members[1:548]
}

# memory.limit(size=12000)   # at least 12GB memory required. 12GB Physical memory strongly recommended

# members matrix
trim <- function (x) gsub("^\\s+|\\s+$", "", x)
make.members.matrix <- function(){
	
	members2 = strsplit(members,split="[.]")
	for (i in 1:length(members2)){
		m = members2[[i]]
		mname = trim(paste(m[5:length(m)],collapse="."))
		members2[[i]] = c(m[1:4],mname)
	}
	t(matrix(unlist(members2),nrow=5))
}
members.mat = make.members.matrix()
colnames(members.mat)=c("type","cong","id","party","name")



######################
## run
######################
##### setup
iter = 200000
thin.n = 200
burn = 100000

# estimation
rc = rollcall(dta, legis.names=members)
save(rc, file="rc.Rdata")

#priors
xp1=rep(0,times = nrow(members.mat))
xp1[members.mat[,"party"]=="D"]=-1
xp1[members.mat[,"party"]=="R"]=1
xpv1 = 1

bp= 0
bpv = 0.04


ptim = proc.time()[3]

results = list()

seed.num = 408393
set.seed (seed.num)
x0 = rep(0, times=nrow(members.mat))
x0[members.mat[,"party"]=="D"] = rnorm(sum(members.mat[,"party"]=="D"),mean=-1, sd=2)
x0[members.mat[,"party"]=="R"]=rnorm(sum(members.mat[,"party"]=="R"),mean=1, sd=2)

prior_setting = list(xp=xp1, xpv=xpv1, bp=bp, bpv=bpv)
	

ideal.est<-ideal(rc, maxiter=iter, thin=thin.n, burnin=burn,,normalize=TRUE, verbose=TRUE, store.item = TRUE,
      priors = prior_setting, startvals=list(x=x0) ) 
	
results[[1]]=ideal.est

# different starting values
#for (i in 1:3){
#	seed.num = 408392 + i
#	set.seed (seed.num)
#	x0 = rep(0, times=nrow(members.mat))
#	x0[members.mat[,"party"]=="D"] = rnorm(sum(members.mat[,"party"]=="D"),mean=-1, sd=2)
#	x0[members.mat[,"party"]=="R"]=rnorm(sum(members.mat[,"party"]=="R"),mean=1, sd=2)
#
#	prior_setting = list(xp=xp1, xpv=xpv1, bp=bp, bpv=bpv)
#	
#
#	ideal.est<-ideal(rc, maxiter=iter, thin=thin.n, burnin=burn,,normalize=TRUE, verbose=TRUE,
#     priors = prior_setting, startvals=list(x=x0) ) 
#	
#	results[[i]]=ideal.est
#	save(results , file = "results.Rdata")
#}
# convergence check in "convergence_check.R" 

save(results , file = "results.Rdata")






tim = proc.time()[3]-ptim
tim/3600

